# coding: utf-8

# 作者：雪山凌狐
# 代码复现：线性回归

"""
Linear Regression: 实现了回归，其中包括线性函数的定义，为什么要用线性函数，loss 的意义，梯度下降的意义，stochastic gradient descent
Use Boston house price dataset.
北京 2020 年房价的数据集，为什么我没有用北京房价的数据集呢？
Boston: room size, subway, highway, crime rate 有一个比较明显的关系，所以就观察关系比较容易
北京的房价：！远近，！房况 ==》 学区！！！！ => 非常贵 特别是海淀区
"""
import random
import numpy as np
import pandas as pd
from sklearn.datasets import load_boston
import seaborn as sns
import matplotlib.pyplot as plt

"""
为什么我没有一行一行解释代码
为什么我没有注释

==> 好的 Python 代码 一定是 "自解释(self description)" 的
"""

dataset = load_boston()
data = dataset['data']
target = dataset['target']
columns = dataset['feature_names']

dataframe = pd.DataFrame(data)
dataframe.columns = columns
dataframe['price'] = target

# 可以看一下生成的数据
# print(dataframe.head())

'''
correlation => 如果一个值的增大，会引起另外一个值一定增大，而且是等比例增大 相关系数就越接近于 1
correlation => 0 就是两者之间没有任何关系
correlation => -1 一个值增大 另外一个值一定减小 而且减小是呈相等比例的
'''
# print(dataframe.corr()) # show the correlation of dataframe variables 展示相关系数
# 还可以查看相关关系的热力图
# sns.heatmap(dataframe.corr())
# plt.show()

# 通过观察相关系数热力图就可以知道哪些特征对于结果的影响最大（正面影响或负面影响）
# 挑出来这些影响较大的特征，形成模型
# 最终取出下面这两个：
# RM：小区平均的卧室个数
# LSTAT: 低收入人群在周围的比例

# 获取我们要用的数据列
rm = dataframe['RM']
lstat = dataframe['LSTAT']


def linear(x, w, b):
    '''
    def model(rm, lstat, w1, w2, b):
        # 我们假设关系是线性的
        return w1 * rm + w2 * lstat + b
    上面的这个函数，就是我们一个线性关系，配上权重和误差的函数。
    我们可以简化成下面正式代码的形式，使用矩阵的点乘的形式来处理。
    我们把 x = (rm , lstat)
          w = (w1, w2)
    这样就可以使用点乘的形式来矩阵相乘，简化模型。

    为什么要抽象成这样的模型呢？
    因为这样写当我们所取的特征变多的时候，点乘依然不用改代码就可以对应位置相乘，具有可扩展性。
    '''
    # vectorized model，w.T 表示向量的转置，转置后才能进行对应位置点乘（参见矩阵的算法知识，比如一行三列跟三行一列才能点乘）
    return np.dot(x, w.T) + b


def loss(yhat, y):
    '''
    最简单的翻译写法：
    def loss(yhat, y):
        loss_ = 0
        for yi, yhat_i in zip(y, yhat):
            loss_ += (yi - yhat_i) ** 2
        return loss_ / len(yhat)
    但是这样写有点代码太多，我们可以用 numpy 提供的求数组均值的方法来实现，达到代码简写的效果。
    '''
    # numpy broadcast，numpy 的广播方法（将一些操作应用于数组中的对应位置的每一个，这里 yhat，y 视为两个数组）
    return np.mean((yhat - y) ** 2)


# 这里是 w1，w2 都给处理了，分别放在元组的两个元素中
def partial_w(x, y, yhat):
    return np.array([2 * np.mean((yhat - y) * x[0]), 2 * np.mean((yhat - y) * x[1])])


# b 的导函数同理
def partial_b(x, y, yhat):
    return 2 * np.mean((yhat - y))


def optimize(w, b, x, y, yhat, pw, pb, learning_rate):
    # 通过学习率获取一个新的 w 和 b，这是梯度下降的方法，获取更优的 w 和 b
    # 有关梯度下降的公式方法，是固定的，可以参见 第二节 机器学习初探 关于梯度下降的解释
    w = w + (-1) * pw(x, y, yhat) * learning_rate
    b = b + (-1) * pb(x, y, yhat) * learning_rate

    return w, b

# 模型训练函数 
def train(model_to_be_train, target, loss, pw, pb):

    # 初始化
    w = np.random.random_sample(size=(1, 2))  # w normal 生成一个一行两列的随机数（正态分布）
    b = np.random.random()  # 0 初始化在 0 附近的随机值 原因深度学习的时候会详细解释
    learning_rate = 1e-5  # 学习率
    epoch = 200  # 时代，指的是在机器学习里面运行了整整一遍
    losses = []

    for i in range(epoch):
        batch_loss = []
        for batch in range(len(rm)):
            # batch training 批训练，实际工作做模型时会随机找几个值来训练模型尝试
            # 因为实际工作中数据量可能非常巨大，一个一个取数验证不是最好实践，随机找会比较好
            # 这种梯度下降的方式叫做 stochastic gradient descent，简称SGD，画出来的图是锯齿状波浪向下的而不是平滑的
            index = random.choice(range(len(rm)))
            rm_x, lstat_x = rm[index], lstat[index]  # 还有一种梯度下降方式是不像这里的每次取一个值来计算，而是一次取多个，叫做 batch gradient descent
            x = np.array([rm_x, lstat_x])
            y = target[index]

            yhat = model_to_be_train(x, w, b)  # 放入模型中
            loss_v = loss(yhat, y)  # 求解 loss 值

            batch_loss.append(loss_v)

            # pw，pb 表示 partial w，partial b
            w, b = optimize(w, b, x, y, yhat, pw, pb, learning_rate)

            # 每训练 100 次做一次显示
            if batch % 100 == 0:
                print('Epoch: {} Batch: {}, loss: {}'.format(i, batch, loss_v))
        # 每一轮求个 loss 平均值
        losses.append(np.mean(batch_loss))

    return model_to_be_train, w, b, losses


if __name__ == "__main__":
    import matplotlib.pyplot as plt

    target = dataframe['price']

    model, w, b, losses = train(linear, target, loss, partial_w, partial_b)
    # 画 losses 的图
    plt.plot(losses)
    # 尝试预测，w b 就是自动求解得到的值了
    predicate = model(np.array([19, 7]), w, b)
    print(predicate)

    plt.show()



